Testing Isopy's three functions using data from Argo float products
import numpy as np
import xarray as xr
from matplotlib import pyplot as plt
%matplotlib inline
url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/Argo_Products/monthly_mean/Gridded_monthly_mean'
ds = xr.open_dataset(url)
ds
<xarray.Dataset>
Dimensions: (lat: 180, lev: 27, lon: 360, time: 180)
Coordinates:
* time (time) datetime64[ns] 2005-01-15 2005-02-15 ... 2019-12-15
* lev (lev) float64 0.0 5.0 10.0 20.0 ... 1.4e+03 1.5e+03 1.75e+03 2e+03
* lat (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
* lon (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
Data variables:
temp (time, lev, lat, lon) float32 ...
salt (time, lev, lat, lon) float32 ...
ptemp (time, lev, lat, lon) float32 ...
pden (time, lev, lat, lon) float32 ...
addep (time, lev, lat, lon) float32 ...
spice (time, lev, lat, lon) float32 ...
Attributes:
title: 1x1 gridded Monthly mean on Standard Levels (from 2005)
Conventions: COARDS\nGrADS
dataType: Grid
documentation: http://apdrc.soest.hawaii.edu/projects/Argo/index.html
history: Fri Jul 19 13:26:16 HST 2019 : imported by GrADS Data Ser...
# %load ../isopy/isopy.py
import numpy as np
import xarray as xr
from scipy import interpolate
def compute_iso(data, iso, lev, d_or_t):
""" compute the depth of some isopycnal/isotherm.
Parameters:
----------------
data: the input variable (Temperature, Density),
array_like (3D), shape (N, M, L).
iso: the objective isosurface layer,
float.
lev: the vertical depth,
array_like (1D), shape (N).
d_or_t: input d or t for density or temperature (lowercase)
string
Returns:
----------------
var_iso: the output variable,
array_like (2D), shape (M, L).
"""
size = data.shape
L = size[2]
M = size[1]
N = size[0] #this should be the same length as lev
#check that the input data is in the for depth lat lon (or depth lon lat)
if len(lev) != N:
print('Error: either data input variable does not have depth as dimension 0, or lev variable is not same length as depth')
else:
if d_or_t in {'d','t'} :
var_iso = np.zeros((M, L)) # define the output var
var_iso[:, :] = np.nan # NaN fill to avoid any later computation errors with zeros
for i in np.arange(L): #loop through dimension 1
for j in np.arange(M): #loop through dimension 2
data_prof = data[:, j, i] #select one profile
if d_or_t == 'd':
id1 = np.where(data_prof < iso)
if np.size(id1) > 0 and np.size(id1) < len(z):
var1 = data_prof[id1[0][-1]]
var2 = data_prof[id1[0][-1] + 1]
if var1 < var2:
func = interpolate.interp1d([var1, var2], [lev[id1[0][-1]], lev[id1[0][-1] + 1]])
var_iso[j, i] = func(iso)
else:
id1 = np.where(data_prof > iso)
if np.size(id1) > 0 and np.size(id1) < len(z):
var1 = data_prof[id1[0][-1]]
var2 = data_prof[id1[0][-1] + 1]
if var1 > var2:
func = interpolate.interp1d([var1, var2], [lev[id1[0][-1]], lev[id1[0][-1] + 1]])
var_iso[j, i] = func(iso)
else:
print('Is this density or temperature?')
return var_iso
def iso_val(density_data, ts_data, iso):
""" compute temperature or salinity along some isopycnal/isotherm.
Parameters:
----------------
density_data: density input,
array_like (3D), shape (N, M, L).
ts_data: temp or salt input,
array_like (3D), shape (N, M, L).
iso: the objective isosurface layer,
float.
Returns:
----------------
ts_iso: the output variable,
array_like (2D), shape (M, L).
"""
size = density_data.shape
L = size[2]
M = size[1]
N = size[0]
ts_iso = np.zeros((M, L)) # define the output var
ts_iso[:, :] = np.nan # NaN fill to avoid any later computation errors with zeros
dens = xr.DataArray(density_data)
var = xr.DataArray(ts_data)
a = dens.where(dens.min(axis=0)<iso)
b = var.where(dens.min(axis=0)<iso)
for j in range(M):
for k in range(L):
x = a[:,j,k]
y = b[:,j,k]
f = interpolate.interp1d(x.values,y.values,'linear')
ts_iso[j,k] = f(iso)
return ts_iso
def iso_average(data, var_iso, lev):
""" compute vertical average above the objective isopycnal layer
Parameters:
----------------
data: the input variable (temperature, salinity, U, V),
array_like (3D), shape (N, M, L).
var_iso: the isopycnal layer,
array_like (2D), shape (M, L).
lev: the vertical depth,
array_like (1D), shape (N).
Returns:
----------------
var_ave: the output variable,
array_like (2D), shape (M, L).
"""
size = data.shape
L = size[2]
M = size[1]
N = size[0]
data_new = (data[:-1, :, :] + data[1:, :, :]) / 2
dz = lev[1:] - lev[:-1]
var_out = np.zeros((M, L))
var_out[:, :] = np.NaN
for i in np.arange(L): #loop through dimension 1
for j in np.arange(M): #loop through dimension 2
var1 = data_new[:, j, i] * dz
idx = find_nearest(lev, dep_iso[j, i])
if idx > 0:
var_out[j, i] = np.nansum(var1[:idx]) / np.sum(dz[:idx])
return var_out
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
compute_iso¶%%time
# read density, temp, salinity
pden = ds.pden #.chunk(chunks={'lat' : 30, 'lon' : 30})
temp = ds.temp
sal = ds.salt
z = pden.lev.values
iso = [27]
for t in np.arange(1):
# isopycnal depth
data0 = pden[90, :, :].values
dep_iso = compute_iso(data0, iso, z, 'd')
# now, try isothermal depth
data1 = temp[90, :, :].values
temp_iso = compute_iso(data1, iso, z, 't')
/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:47: RuntimeWarning: invalid value encountered in less /srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:56: RuntimeWarning: invalid value encountered in greater
CPU times: user 5.18 s, sys: 173 ms, total: 5.35 s Wall time: 13 s
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def make_map(projection=ccrs.PlateCarree(), figsize=(20, 15)):
fig, ax = plt.subplots(
figsize=figsize,
subplot_kw={"projection": projection})
return fig, ax
fig, ax = make_map(projection=ccrs.PlateCarree(), figsize=(20,15))
ax.set_global()
ax.coastlines(resolution="110m", color="k")
ax.add_feature(cfeature.LAND, facecolor="0.75")
lat = ds.lat.values
lon = ds.lon.values
lons, lats = np.meshgrid(lon, lat)
h = ax.contourf(lons, lats, dep_iso, trans=ccrs.PlateCarree())
ax.contour(lons lats, dep_iso, trans=ccrs.PlateCarree(), colors='w')
plt.colorbar(h, orientation='vertical', fraction=0.02).set_label(label='ËšC', size=15)
ax.set_title('Depth field of the \nabove 90 isopycnal mark', size=20);
/srv/conda/envs/notebook/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1400: UserWarning: The following kwargs were not used by contour: 'trans' result = matplotlib.axes.Axes.contourf(self, *args, **kwargs) /srv/conda/envs/notebook/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: The following kwargs were not used by contour: 'trans' result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
# change numpy to Xarray
ds_dep_iso = xr.DataArray(dep_iso, coords={'lat': lat, 'lon': lon}, dims=['lat', 'lon'])
ds_dep_iso
<xarray.DataArray (lat: 180, lon: 360)>
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]])
Coordinates:
* lat (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
* lon (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
import hvplot.xarray
import geoviews as gv
import holoviews as hv
from holoviews import opts
proj = crs.PlateCarree()
ds_dep_iso.hvplot.quadmesh(
'lon', 'lat', projection=proj, project=True, global_extent=True,
width=1000, height=540, cmap='magma_r', rasterize=True, dynamic=False
) * gv.feature.coastline * gv.feature.land(color='red')
# ds_dep_iso.hvplot.contour('lon', 'lat', projection=proj, project=True, global_extent=True,
# width=1000, height=540, rasterize=True, dynamic=False, line_width=5
# ) * gv.feature.coastline
iso_val¶%%time
isotemp = iso_val(pden[90, :, :].values, temp[90, :, :].values, iso)
isosal = iso_val(pden[90, :, :].values, sal[90, :, :].values, iso)
/srv/conda/envs/notebook/lib/python3.7/site-packages/xarray/core/nputils.py:220: RuntimeWarning: All-NaN slice encountered result = getattr(npmodule, name)(values, axis=axis, **kwargs)
fig, ax = make_map(projection=ccrs.PlateCarree(), figsize=(20,15))
ax.set_global()
ax.coastlines(resolution="110m", color="k")
ax.add_feature(cfeature.LAND, facecolor="0.75")
lat = ds.lat.values
lon = ds.lon.values
lons, lats = np.meshgrid(lon, lat)
h1 = ax.contourf(lons, lats, isotemp, np.arange(0,16,0.5), trans=ccrs.PlateCarree(), cmap='RdBu',
vmin=2, vmax=14)
plt.colorbar(h1, orientation = 'vertical', fraction=0.02, label='ËšC').set_label(label='ËšC',
#weight='bold',
size=15)
ax.set_title('Temperature field along the \n 90 isopycnal mark', size=20);
iso_average¶temp_avr = iso_average(temp[90, :, :].values, temp_iso, ds.lev.values)
temp_avr
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]])
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.axes_grid1 import AxesGrid
def make_map(projection=ccrs.PlateCarree(), figsize=(20, 15)):
fig, ax = plt.subplots(
figsize=figsize,
subplot_kw={"projection": projection})
return fig, ax
fig, ax = make_map(projection=ccrs.PlateCarree(), figsize=(20,15))
# fig, ax = plt.plot(projection = ccrs.PlateCarree(), figsize = (20,15))
# fig = plt.figure(figsize=(20,15))
# ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines(resolution="110m", color="k")
ax.add_feature(cfeature.LAND, facecolor="0.75")
lat = ds.lat.values
lon = ds.lon.values
#lons, lats = np.meshgrid(lon, lat)
h2 = ax.contourf(lon, lat, temp_avr, np.arange(0,20,2.0), extend='both', trans=ccrs.PlateCarree(),
cmap='RdBu')
plt.colorbar(h2, orientation = 'vertical', fraction=0.02).set_label(label='ËšC',
#weight='bold',
size=15)
ax.set_xticks(np.linspace(-180, 180, 5), ccrs.PlateCarree())
ax.set_yticks(np.linspace(-90, 90, 5), ccrs.PlateCarree())
ax.set_title('Vertically averaged temperature\nabove the 90 isopycnal mark', size=20);
# lon_formatter = LongitudeFormatter(zero_direction_label=True)
# lat_formatter = LatitudeFormatter()
# ax.xaxis.set_major_formatter(lon_formatter)
# ax.yaxis.set_major_formatter(lat_formatter)